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Abstract 

Background: Ihe expression of myogenic regulatory factors (MRFs) consisting of MyoD, Myf5, myogenin {MyoC) and MRF4 
characterizes various phases of skeletal muscle development including myoblast proliferation, cell-cycle exit, cell fusion and 
the maturation of myotubes to form myofibers. Although it is well known that the function of MyoG cannot be 
compensated for other MRFs, the molecular mechanism by which MyoC controls muscle cell differentiation is still unclear. 
Therefore, in this study, RNA-Seq technology was applied to profile changes in gene expression in response to MyoC knock- 
down (IVlyoGkd) in primary bovine muscle satellite cells (MSCs). 

Resu/ts: About 61-64% of the reads of over 42 million total reads were mapped to more than 1 3,000 genes in the reference 
bovine genome. RNA-Seq analysis identified 8,469 unique genes that were differentially expressed in MyoGi^j. Among these 
genes, 230 were up-regulated and 224 were down-regulated by at least four-fold. DAVID Functional Annotation Cluster 
(FAC) and pathway analysis of all up- and down-regulated genes identified overrepresentation for cell cycle and division, 
DNA replication, mitosis, organelle lumen, nucleoplasm and cytosol, phosphate metabolic process, phosphoprotein 
phosphatase activity, cytoskeleton and cell morphogenesis, signifying the functional implication of these processes and 
pathways during skeletal muscle development. The RNA-Seq data was validated by real time RT-PCR analysis for eight out of 
ten genes as well as five marker genes investigated. 

Conclusions: Jh'\s study is the first RNA-Seq based gene expression analysis of MyoG^d undertaken in primary bovine MSCs. 
Computational analysis of the differentially expressed genes has identified the significance of genes such as SAP30-like 
(SAP30L), Protein lyl-1 (LYLl), various matrix metalloproteinases, and several glycogenes in myogenesis. The results of the 
present study widen our knowledge of the molecular basis of skeletal muscle development and reveal the vital regulatory 
role of MyoC in retaining muscle cell differentiation. 
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Introduction 

Skeletal muscle formation is a multi-step process that requires 
proliferation of myocytes, expression of muscle-specific myogenic 
regulatory factors (MRFs) including MyoD, MyJ5, myogenin [MyoG] 
and MRF4 (or Myf6), cell cycle withdrawal, myotube formation by 
the fusion of mononucleated cells and maturation of myotubes into 
myofibers [1], [2], [3], [4], [5]. MRFs are basic helix-loop-helix 
(bHLH) transcription factors [6] that cooperate with several 
transcription factors of the myocytes enhancer factor-2 {MEF2) family 
[7] to regulate myogenesis. bHLH proteins also heterodimerize 
with E-proteins [8], enabling binding to the E-Box consensus 
sequence (CANNTG) present in the regulatory regions of muscle 



specific genes [9], [10]. Among these MRFs, MyoD is highly 
expressed during the mid-Gl phase and between the S and M 
phases of the cell cycle, but absent during the GO phase [11], 
whereas Myf5 is highly expressed during the GO phase and 
decreases during the Gl phase [12]. MyoG and MRF4 {Myf6) are 
expressed upon differentiation of myoblasts into multinucleated 
myotubes [13], [14], [15]. 

MyoG is crucial during differentiation [11], as many studies have 
revealed that mice lacking MyoG continue to identify the muscle 
lineage through the formation of myoblasts [16], but show high 
perinatal mortality due to severe skeletal muscle deficiency caused 
by disruption of myoblast differentiation and muscle fiber 
formation [17], [18]. Additionally, MyoG/ MyoD and MyoG/MyP 
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double knockout mice studies have shown that these mice specify 
the muscle lineage, but the formation of muscle fibers is disrupted, 
which is similar to MyoG knockout mici; [19]. Furthermore, MyoD 
and Myf5 are unable to compensate for the role of AdyoG in 
differentiation [20], and mice that lack AdyoG exhibit normal 
expression levels oi MyoD and A'lyf5 [17]. This is because MyoG 
acts downstream of MyoD and Myf5 [16] in skeletal muscle 
differentiation. Knockout mice studies have also shown a 
relationship between different MRFs in which the absence of 
one will be compensated for by another [21], [22]. The only 
exception to this compensation effect of MRFs is MyoG, which 
plays a unique and non-redundant role during embryogenesis 
[19], whereas conditional knock-out resulted in reduced muscle 
mass in adults [23]. 

In this study, we conducted a comprehensive transcriptome 
analysis of primary- bovine cells using MyoG^ and compared the 
expression profiles with those of the wild type using an RNA-Seq 
tedmique. We also showed that MyoGkd led to upregulation of 
genes involved in processes such as cell proliferation and DNA 
replication, whereas the genes involved in phosphate metabolic 
processes were down-regulated. Finally, potential involvement of 
various new genes in myogenesis was identified. 

Materials and Methods 

Bovine MSCs culture 

Bovine muscle was collected from the hind leg skeletal muscle of 

24-26 month old cattle with a body weight of 550-600 kg. The 
animals were handled according to a protocol approved by the 
Animal Care and Concern Committee of the National Institute of 
Animal Science, Korea. Briefly, the collected muscle was minced 
into fine pieces, and digested with trypsin-EDTA (GIBCO, CA, 
USA), and were centrifuged at 90 xg for 3 min and the upper 
phase was passed through a 40-nm cell strainer. The filtrate was 
centrifuged at 2,500 rpm, pellet was collected, washed twice and 
cultured in Dulbecco's modified Eagle's medium (DMEM; 
HyClone Laboratories, UT, USA) supplemented with lO'M) fetal 
bovine serum (FBS, HyClone Laboratories) and 1% penicillin/ 
streptomycin at 37°C under 5% COg. The culture medium was 
changed every other day. To induce differentiation, cells were 
allowed to grow in DMEM without reducing serum (DMEM with 
10% FBS and 1% P/S) for 10, 12, 14, 16, and 18 days. MSCs 
isolation and culture were conducted as previously described [24] . 

MyoG shRNA construction and knock-down 

Bovine MyoG shRNA was designed using nucleotide informa- 
tion obtained from NCBI (AB257560.1) and cloned with pRNAT- 
U6.2/Lenti vector (GeneScript, NJ, USA). Constructed MyoG 
shRNA or non-specific sequences (scrambled vector, MyoG„d) 
were transfected to generate viral particles in 293 FT cells. After 
two days of transfection, the supernatant containing viral particles 
was collected, transduced with lentiviral particles expressing 
shRNAs against bovine MyoG or scrambled vector in MSCs 
(Day 8), and selected with 50 ng/ml of G418 (CABIOCHEM, 
CA, USA). The selected cells were allowed to differentiate and 
were harvested at Day 2 1 . The following oligonucleotide was used 
to generate AdyoG %hKNA: sense: 5'- GGATCCCGCGCAGACT- 
CAAGCCGCCGGTGTTCAAGAGACACCTTCTTGAGTCTG- 
CGCTTTTCCAACTCHGAG-3'. 

RNA extraction, library preparation and sequencing 

MSCs were allowed to grow till day 10, and were transduced 
with either scrambled vector or MyoG shRNA. Cells were then 
allowed to grow for another 1 1 days, and were harvested with 



Trizol reagent (Invitrogen) according to the manufacturer's 
protocol. Total RNA was then extracted and stored in 
diethylpyrocarbonate-treated HgO at — 80°C until used. The 
mRNA in total RNA was converted into a library of template 
molecules suitable for subsequent cluster generation using the 
reagents provided in the lUumina TruSeq RNA Sample Prepa- 
ration Kit (lUumina, CA, USA) according to the manufacturer's 
instructions. Library construction and high-throughput sequencing 
were carried out using an lUumina HiSeq2000 sequencing system 
in which each sequencing cycle takes place in the presence of all 
four nucleotides, leading to higher precision than methods in 
which only a single nucleotide is present in the reaction mixture at 
one time. The cycle is repeated one base at a time, creating a 
string of images each indicating a single base extension at a specific 
cluster. 

Sequence quality check 

The FASTQC [http://www.bioinformatics.babraham.ac.uk/ 
projects/fastqc/] tool embedded in the web-based platform. 
Galaxy [25], [26], [27], was used to calculate quality control 
statistics describing raw sequence data from FASTQ^ files 
generated by the lUumina second generation sequencing technol- 
ogy ("Solexa") [http://www.illumina.com/technology/solexa_ 
technology.iknn] . 

Mapping of RNA-Seq reads transcript assembly 

TopHat [28] was used to ahgn RNA-Seq reads against UCSC 
Bos taums reference genome (Btau_4.6.1/bosTau7) via Bowtie, 
which is a very high-throughput short read aligner [29] . Bowtie is 
different from other general-purpose alignment tools such as 
BLAST [30], and shows best performance when short reads are 
ahgned to large genomes. Bowtie is extremely fast for short reads 
where several reads have at least one significantly valid alignment, 
the reads are of high-quality, and the number of alignments 
reported pf;r r(;ad is nearly 1 [29], These mapping results were 
then analyzed to identify splice junctions between exons. AH 
default parameters were used to run TopHat except the mate 
inner distance, for which a value of 100 was selected in the case of 
paired reads. The advantage of a paired end run is that both reads 
contain long range positional information, allowing for highly 
precise alignment of reads. 

The aligni;d reads were further analyzed by Cufflinks [31] using 
a multifasta file (bosTau7. fa) option that can improve the 
precision of transcript abundance approximation by bias detection 
and a correction algorithm. The relative abundance of transcripts 
was reported as fragments per kilobase of exon per miUion 
fragments mapped (FPKM). An additional cufflinks parameter for 
the initial estimation procedure was used so that the reads 
mapping to multiple locations in the genome were accurately 
weighted [31]. The nucleotide sequences obtained in this study 
have been submitted and will be available in NCBI Short Read 
Archi\T with accession number SRR 1122446 as soon as it is 
released. Alternatively, the data can be obtained directly from the 
authors. 

Functional annotation cluster and pathway analysis 

DAVID [http:/ /david.abcc.ncifcrf g()\/h<)me.jsp] functional 
annotation cluster analysis was performed on the hst of up- 
regulated and down-regulated genes with a fold change of &4. 
Only those terms that reported a jd-value of £0.05 and count 
number a 5 genes were selected for analysis. The Gene Ontology 
(GO) terms of cellular component, molecular function and 
biological process in DAVID were employed to categorize 
enriched biological themes in up- and down-regulated gene lists. 
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Pathway mapping was performed using the KEGG Automatic 
Annotation Server (KAAS) [32]. The nucleotide sequences of up- 
and down-regulated gc-ncs were uploaded to the KAAS web server 
as an input using single-directional best hit (SBH) metliod to assign 
orthologs. KAAS offers functional annotation of genes in a 
genome via a BLAST similarity searches against a manually 
curated set of ortholog groups in the KEGG GENES database. 
KAAS assigned a KEGG Orthology (KO) number to genes in the 
data sets, which were mapped to one of KEGG's reference 
pathways. 

Real time RT-PCR validation 

One microgram of RNA in a reaction mixture with a total 
volume of 20 |^1 was primed with oligo (dT)2() primers (Bioneer, 
Daejeon, Korea) and then reverse transcribed at 42°C for 50 min 
and 72°C for 15 min. Subsequently, 2 |il of cDNA product and 10 
pmoles of each gene-specific primer were used for PGR, using a 
7500 real-time PGR system (Apphed Biosystems, Foster City, CA, 
USA). A Power SYBRH Green PGR Master Mix (Applied 
Biosystems) was used as the fluorescence source. Primers were 
designed with the Primer 3 software (http://frodo.wi.mit.edu) 
using the sequence information hsted at the National Center for 
Biotechnology Information. Detailed information describing the 
primer sequences is provided in Table SI. 

Immunocytochemistry 

Cells grown in a covered glass-bottom dish were stained with 
Pax7 or MyoG antibody. Briefly, cells were rinsed with PBS 
(phosphate buffered saline), fixed in 4% formaldehyde, permea- 
bilized by 0.2% TritonX-100, after which the signals were 
enhanced using an Image-iT FX signal enhancer (Invitrogen). 
The cells were then incubated with mouse primary Pax7 or MyoG 
antibody (1 :50, Santa Cruz Biotechnology, CA, USA) at 4°C in a 
humid environment overnight. Secondary antibody (Alexa Fluor 
488 goat anti-mouse SFX kit; Molei:ular Probes, Eugene, OR, 
USA) was treated for 1 hr at room temperature followed by 
nuclear staining with 4',6'-diamino-2-phenylindole (DAPI; Sigma- 
Aldrich, MO, USA). Pictures were taken using a fluorescent 
microscope equipped with a digital camera (Nikon, Tokyo, Japan). 

Western blot 

Western blot was performed with the total protein isolated from 
cells. Briefly, cells washed with ice-cold PBS were lysed in RIPA 
lysis buffer containing protease inhibitor cocktail (Thermo 
Scientific, IL, USA). The protein was quantified by Bradford 
method using protein assay dye solution. Fifty microgram of 
protein was electrophoresed in 10% SDS-polyacrylamide gel after 
reducing at 90°C for 3 min with P-mercaptoethanol, and the 
protein was transferred to a PVDF membrane. Membrane was 
blocked and hybridized with MyoG (1:1000) or P-actin antibody 
(1:2000) (Santa Cruz Biotechnology) ()\ <;rnight at 4' C. Membrane 
washed in TBST was then incubated with horseradish peroxidase 
conjugated secondary antibody for an hour at room temperature. 
Finally, the membrane was developed using SuperSignal West 
Pico Chemiluminescent Substrate (Thermo Scientific). 

Giemsa staining 

Cells were washed with PBS, fixed with PBS/methanol (v/v) for 
2 min, and were incubated with 0.04% Giemsa G250 solution for 
30 min. cells were rinsed with distilled water and pictures were 
taken using a light microscope equipped with a digital camera 
(Nikon). 



Statistical analysis 

The normalized expression means were compared using 
Tukey's Studentized Range (HSD) to identify significant differ- 
ences in gene expression. A nominal /)-value of less than 0.05 was 
considered to be statistically significant. Real time RT-PCR data 
were analyzed by one-way ANOVA using PROG GLM in SAS 
package ver. 9.0 (SAS Institute, Gary, NC, USA). 

Results 

MyoG gene knock-down 

MSCs isolated and cultured from bovine leg muscle were 
stained with Pax7 to determine the purity of cells. The isolated 
cells showed approximately 85% Pax7 positive cells (Figure SI). 
The in vitro cultured bovine primary cells began to difierentiate 
without serum deprivation. The initial myotubes became visible at 
Day 12, and the number of myotubes increased with time 
(Figure lA). The expression level oiMyoG, which is known to play 
a role in myogenic differentiation [33], was determined by real- 
time RT-PCR at different time points of primary bovine cells 
difierentiation. MyoG was expressed throughout difierentiation; 
however, there was a steep increase in its expression level at Day 
12, which gradually increased until Day 18 (Figure IB). 
Similarly, Western blot analysis revealed MyoG protein expression 
on Day 10 and Day 16 with higher levels occurring on Day 16 
(Figure IC). This expression profile of MyoG during differenti- 
ation is in accordance with those of previous studies [34], [35]. 
Moreover, the cells were authenticated to be in the state of 
myotube formation by inspecting the nuclear expression of MyoG 
protein at two different time points (Day 10 and Day 16) during 
cell differentiation. MyoG protein expression was observed at Day 
10 and 16. Day 16 showed higher MyoG nuclear expression as 
compared with Day 10 proliferating cells (Figure ID). To identify 
the genes differentially expressed as a consequence of MyoG knock- 
down, MSCs were transduced with sliRNA spc-cific for MyoG. 
RNA analysis following transduction rc\ (-aled the specific decline 
of mRNA for shRNA induced MyoG^i as compared to its wild 
type counterpart (MyoG^a) (Figure IE). Similarly, MyoGkd was 
confirmed at the protein level by Western blot analysis 
(Figure IF). shRNA transduction against MyoG prohibited the 
nuclear expression of MyoG protein and the development of 
myotubes (Figure IG and IH). 

Expression of myogenic marker genes 

To confirm that primary bovine cells were undergoing 
differentiation, we verified the expression of myosin regulatory light 
chain 2 {MYL2] and myosin heavy chain 3 {MTH3}, which have 
previously been shown to be expressed during myogenesis. Both 
MTL2 and MYH3, which are marker genes [24], exhibited a 
gradual increase in expression rates during myogenesis, whereas 
cyclin A2 [CCNA2), which is involved in the cell cycle [36], showed 
moderate and decreased expression levels (Figure 2A). The 
opposite trend was observed for MYL2 and MYH3, with decreased 
mRNA expression, while the expression of CCMA2 was signifi- 
cantiy elevated as a result of MyoGjsd (Figure 2B). These results 
are in agreement with those of our previous study [37], [38], as 
well as those of other investigations of mouse and human skeletal 
muscle differentiation [4], [39]. 

High-throughput sequencing 

High-throughput RNA-Seq was applied to investigate the gene 
expression profiles of MyoG„,| and MyoGj;,! samples. The total 
numbers of RNA-Seq reads (101 base pairs in length) generated in 
this study were about 42 million for MyoG„d and 46 million for 
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Figure 1. Myogenesis is associated witKi increased MyoG expression. A) MSCs proliferation and differentiation in DiVIEIVI/10% FBS media. 
IVIyotube formation was observed from day 1 2, reaching a maximum at day 1 6. B) MyoG expression from RNA extracted and analyzed by real time RT- 
PCR. MyoG expression gradually increased with time. C) IVlyoG protein expression at Day 1 0 and Day 1 6 was determined by Western blot analysis and 
the protein intensity was measured by using Image J program. IVlyoG protein in Day 16 cells was significantly higher than in Day 10 cells. D) 
Immunocytochemistry analysis of Day 1 0 and Day 1 6 cells stained with MyoG antibody. Day 1 6 showed higher IVlyoG nuclear expression as compared 
with Day 10 proliferating cells. E)& F) MyoG mRNA and protein expression of MSCs transfected with MyoG shRNA. MyoG expression was decreased in 
both RNA and protein level in MyoG^d cells relative to MyoG„d cells. G) Immunostaining of MyoG^j or MyoG^d cells with MyoG antibody. A significant 
decrease in nuclear MyoG protein expression was observed in MyoG|<d cells. H) Giemsa staining performed in MyoG„;j and MyoG^d cells. Myotube 
formation was decreased in MyoGkd cells relative to MyoG„d cells. Day 10 and MyoG„d represents control, respectively (mean ±S.D., n = 3). p-value 
indicates the statistical significance of the data and different letters (a and b) in graph show significant differences among groups. 
doi:1 0.1 371 /journal.pone.0092447.g001 



MyoGkd. About 63.73% and 61.66% of the MyoG„d and MyoGkd 
reads, with at least one reported alignment, were mapped to the 
reference genome (Table 1). Post-run quality analysis of the 
RNA-Seq data was carried out using the FASTQG [http:/ /www. 
bioinformatics.babraham.ac.uk/projects/fastqc/] tool in Galaxy 
[25], [26], [27]. The per base sequence quality report is one of the 
most useful FASTQC reports, which helps in deciding whether 
sequence trimming is required before alignment. Figure 3A— D 
summarizes the range of quality standards across all bases at every 
point in the FastQ^fde. For each position a boxwhisker type plot is 
drawn. In general, the quality of calls will degrade as the run 
advances, therefore, it is prevalent to see base calls falling into the 
orange area towards the end of a read. The quality scores across 
all bases were determined by the Sanger/IUumina 1.9 encoding 
method. These figures represent good quality calls scattered across 
the green background of the plots. A warning wiU be pointed out if 
the lower quartUe for any base is less than 1 0, or if the median for 
any base is less than 25, whereas a failure will be issued if the lower 
quartUe for any base is less than 5 or if the median for any base is 
less than 20. 



Differentially expressed genes 

Following mapping of the sequencing reads to the reference 
genome with TopHat [28], transcripts were assembled and their 
relative expression levels were computed with Cufflinks [31] in 
FPKM. A total of 13,703 unique genes were detected and further 
filtered to remove possible noise from the data by excluding the 
genes with FPKM values equal to zero from the analysis. As a 
result, 9,337 and 12,835 genes were identified from MyoG„d and 
MyoGi^d samples respectively, which shared 8,469 genes in 
common (Table 2). These 8,469 genes were then used to 
calculate the fold change, which was defined as the ratio of 
MyoGkd FPKM to MyoG„d FPKM. In this study, the total fold 
change of S4 was considered to classify the differentially expressed 
genes. Based on this defmition, there are 230 up-regulated and 224 
down-regulated genes in MyoG^d over the MyoG„d sample 
(Table S2). We found that SAP30-like {SAP30L) was the most up- 
regulated gene in MyoGkd by 126-fold (logs fold change =6.98) 
and encodes a protein that plays a potential role in the histone 
deacetylase complex, similar to Sin3 associated protein 30 (SAP30) 
[40]. 
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Figure 2. Effect of MyoGkd on MYL2, MYH3 and CC/V/12 genes. A) Time course study of mRNA expression of MYL2, MYH3 and CCNA2 during 
IVlSCs differentiation. Cells cultured in differentiation media showed gradual increase in MYL2 and MYH3 expression until Day 18, but transient 
increase of CCNA2 at day 1 2 decreased at later stage of myogenesis. B) Evaluation of MYL2, MYH3, and CCNA2 gene expression by real time RT-PCR in 
MyoGkd cells. Decreased MYL2 and MYH3 gene expression and increased CCNA2 expression was observed in IVlyoG|<d cells relative to MyoG|<d cells. Day 
10 and MyoG„d represents control, respectively (mean ±S.D., n = 3). p-value indicates the statistical significance of the data and different letters (a, b 
and c) in graph show significant differences among groups. 
doi:1 0.1 371 /journal.pone.0092447.g002 



SAP30 was identified as one of the transcriptional regulators in 
C2CI2 difiFerentiation [41]. Ribosomal protein L23a {RPL23A), zinc 
finger protein 322A (ZNF322A), solute carrier fiimily 16 member 3 
{SLC16A3), tubulin, alpha Ic [TUBAIC), sulfotransferase fiimily, cytosolic, 
2B, member 1 (SULT2B1), metallothionein 2A {MT2A), matrix 
metallopeptidase 9 (MMP9), secreted finzzled-related protein 1 (SFRPl) 
and myelin protein (MBP) were among the ten most up-regulated 
genes in MyoG^d. MT2A, a member of cysteine-rich and metal 
binding intracellular proteins [42] that has been linked with cell 
proliferation [43], was up-regulated by 17.77-foId. In addition to 
MT2A, the other two metaUothioneins iMTlA and MTIE] present 
in the list of differentially expressed genes also showed up- 
regulation by more than four-fold. Another group of highly up- 
regulated genes belonged to a class of matrix metalloproteinases 
(MMPs). A total of 13 MMPs were identified in this study, almost 
all of which showed a high fold change. MMPl, MMP3, MMP9 
(one of the ten most up-regulated genes) and MMPl 3 were up- 
regulated by more than four-fold, whereas MMP2 and MMP12 
were up-regulated by more than two-fold. Similarly, MMP15, 
MMP19, MMP20, MMP23B and MMP27 were up-regulated by 
more than one-fold. Only one of the MMPs, MMPl 6, was down- 
regulated. Increased expression of MMPs (particularly MMP9) in 
skeletal muscles is well known [44]; therefore, our analysis is 
consistent with the results of previous studies. These MMPs enable 
release of the active hepatocyte growth factor [HGF), which stimulates 
proliferation while inhibiting differentiation, from extracellular 
matrix (ECM) [45]. The other genes that showed a greater than 
four-fold increase in expression include a transcription repressor, 
musculin also known as MyoR (myogenic repressor), which is known 
to block myogenesis and the activation of E-box dependent muscle 
genes [46]. 

Protein lyl-1 (LTLl), also known as lymphoblastic leukemia-derived 
sequence 1, was found to be one of the ten most down-regulated 
genes. LYLl consists of a basic helix-loop-helix (bHLH) domain 



(Figure S2), which is similar to genes involved in mammalian 
myogenesis [MyoD, MyoG, MyfiS, and herculin) [47]. LTLl is an 
essential gene required for the development of adult hematopoietic 
stem cells [48] . The other ten most down-regulated genes include 
sodium channel protein type 1 subunit alpha {SCMIA), ribosomal protein 
SI 5a {RPS15A), syncoilin {STJVQ, tubulin alpha- ID {TUBA W), family 
with sequence similarity 65, member B (FAM65B], agouti- signaling protein 
[ASIP], tocopherol (alpha) transfer protein-like (TTPAL), ryanodine receptor 1 
[RYRl), and an uncharacterized protein (LOG 100847946). Myocyte 
enhancer factor 2C [MEF2C) is a member of the MEF2 family of 
transcription factors that was down-regtilated by two-fold, whereas 
another member of the MEF2 family, MEF2A, was down- 
regulated by about 1.5-fold. MEF2 transcription factors play 
prominent roles in skeletal muscle differentiation [49], [50], and 
four MEF2 isoforms {MEF2A, MEF2B, MEF2C and MEF2D) have 
been identified to date [51], [52]. It is well known that increased 
expression of occurs during myoblast differentiation [53], 

[54], [55]. In comparison with the up-regulated genes, down- 
regulated genes consisted of those that play crucial roles in 
phosphate metabolic processes. There is a significant amount of 
evidence that the processes related to phosphorylation and 
dephosphorylation of tyrosine are important regulatory compo- 
nents during the progression of myogenesis [56], [57], [58], [59]. 

Functional annotation cluster and pathway analysis 

To categorize biological processes that are overrepresented in 
MyoG„d and MyoGkd cells, we classified all known difTerentially 
expressed genes (fold change &4) using the Functional Annotation 
Gluster (FAC) tool available in the Database for Annotation, 
Visualization and Integrated Discovery (DAVID) [http://david. 
abcc.ncifcrf.gov/home.jsp]. DAVID FAC analysis of 230 up- 
regulated genes (fold change &4) generated a total of 65 functional 
clusters using default parameters. The GO terms "Biological 
Process", "Cellular Component" and "Molecular Function" were 
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Table 1. Number of single replicate RNA-Seq reads across MyoG^j and MyoG^d samples. 





Sample 


Total read pairs° 


Processed reads'* 


Mapped reads'^ 


WlyoGwd 


42,936,654 


42,910,546 


27,346,293 (63.73%) 


MyoGkd 


46,349,131 


46,324,224 


28,565,046 (61.66%) 



Foot note: a. Total read count b. Reads used in TopHat process c. Reads with at least one reported alignment. 
doi:l 0.1 371/journal.pone.0092447.t001 



used for annotations. Genes with a variety of GO terms from tlie 
resulting functional clusters having statistically significant p-vahies 
are listed in Table 3A. Similarly, the GO functional annotation 
chart reported 50 chart records that were further filtered to 22 
records by selecting only those terms having ^-values ^0.05 and 
number of genes in each term &5 (Table S3A). From these tables, 
it can be seen that the GO terms are enriched in genes with 
functions necessary for actively proliferating cells such as cell 
division, DNA replication, cell cycle function and mitosis. The 
other major processes that exhibit higher levels of gene expression 
as a consequence of MyoGj^j include GO terms related to 
organelle lumen, nucleoplasm and cytosol. These data suggest that 
the proliferating processes are the major processes up-regulated by 
MyoGj^a, and that their overrepresentation may be due to the de- 
differentiation of muscle cells. The down-regulation of MyoG in 
terminally differentiated mouse C2C12 myotubes was recendy 



shown to stimulate cellular cleavage into mononucleated cells and 
promote cell cycle re-entry [1 1]. This phenomenon of dediflFeren- 
tiation of myotubes into proliferating mononucleated cells is well 
known [60], [61], [62]. 

Functional analysis of 224 down-regulated genes resulted in 69 
clusters, and the statistically significant (/(-value £0.05) GO terms 
having at least five members in each enriched term are listed in 
Table 3B. The GO functional annotation chart reported 59 
records, 21 of which were selected on the basis of a /(-value £0.05 
and number of genes in each term a5 (Table S3B). The 
processes that were significandy down-regulated by MyoG];;) as 
shown by functional analysis include phosphate metabolic 
processes, dephosphorylatioii, phosphoproteiii phosphatase activ- 
ity and protein amino acid dephosphorylatioii. Additionally, 
processes related to the cell shape such as cytoskeleton and cell 
morphogenesis (Table 3B & Table S3B) were also down- 




Figure 3. Per base sequence quality of iVlyoG„d and IVlyoG^d. Quality scores of A) & B) MyoG^j, and C) & D) IVlyoG|<d. The y-axis on the graph 
shows the quality scores with higher scores indicating better base calls. The background of the graph separates the y axis into high-quality calls 
(green), calls of reasonable quality (orange), and calls of poor quality (red). In each of these findings, the red line is the median value, the yellow box 
corresponds to the inter-quartile range (25-75%), the upper and lower whiskers represent the 10% and 90% points, respectively, and the blue line 
signifies the mean quality. 
doi:1 0.1 371/journal.pone.0092447.g003 
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Table 2. Gene expression summary. 




Sample MyoG„d 


MyoGkd 


Total No. of genes 9,337 


12,835 


Common in MyoG^d/'V'yoGkd 


8,469 


Up-regulated (2:4 fold) 


230 


Down-regulated (a4 fold) 


224 



doi:1 0.1 371 /journal.pone.0092447.t002 



regulated. In addition to DAVID FAC, by performing KEGG 
pathway analysis of 455 up- and down-regulated genes, we were 
able to assign 281 unique KEGG orthologs to these dilferentially 
expressed query genes. The majority of the dififerentiaUy expressed 
genes were found to be associated with important biological 
processes, many being classified in signaling pathway or being 
involved in adhesion and cytoskeleton related functions (Figure 
S3A-S3E). 

Glycogenes expression 

To further explore the role of glycogenes in myogenesis, all 230 
up- and 224 down-regulated genes were manually verified in the 
UniProt database [63] to check whether they represent a 
glycogene or not. If a gene encoded a protein and represented 
glycosyltransferases, glycosidases, lectins, sulfotransferases or 
proteins involved in carbohydrate metabolism or transport [5], it 
was labeled as a glycogene. In this way, we identified 59 (~25%) 
up- and 52 (~23%) down-regulated glycogenes out of 230 and 224 
differentially expressed gene sets. Some of the glycogenes that 
demonstrated four-fold increase in expression rates included 
SCMNIA, SFRPl and transmembrane protein 217 {TMEM217). 
Additionally, various glycogenes such as matrix metalloproteinases 
{MMPl, MMP13, MMP3 and MMP9) and genes belonging to the 
solute carrier family {SLC26A8, SLC2A6, SLC37A2 and SLC46A3) 
exhibited more than four-fold up-regulation as a consequence of 
MyoGjjj. Similarly, the glycogenes that showed four-fold decrease 
included SCNIA, ASIP and protein wnt-11 {WNTll). The hst of top 
ten up- and down-regulated genes consisted of at least two 
glycogenes. 

Validation of RNA sequencing data 

To validate the RNA-Seq results, we performed real-time RT- 
PCR to determine the expression levels of marker genes {Myf5, 
MjioD, MyoG, MTL2 and MYH3) involved in myogenesis and then 
compared their expression with RNA-Seq data in MyoGj^^j 
samples. The RT-PCR results were well correlated with the 
RNA-Seq expression data for the five marker genes investigated 
(Figure 4A). SpecificaUy, RT-PCR analysis oiMjfS, MjoD, MyoG, 
MrL2 and MrH3 mRNA levels revealed fold changes of 0.72, 
1.32, 0.42, 0.35 and 0.67, respectively (approximately >2-fold 
decrease) in MyoGj^a relative to MyoG^j cells. These results 
compliment favorably well with our RNA-Seq data showing fold 
changes of 0.8 {Mjf5), 2.23 {MjoD), 0.49 {MyoG), 0.33 {MrL2) and 
0.49 {MTH3), which also corresponded to a two-fold decline in the 
expression of these marker genes. As an additional confirmation of 
the expression data, ten genes were randomly selected for RT- 
PCR analysis, five of them representing the ten most-up and 
down-regulated genes. The results revealed that the fold-change 
profiles measured by RNA-Seq and RT-PCR were concordant for 
all ten genes. However, RT-PCR analysis showed statistically 
significant expression of only eight out of these additional ten 
genes (Figure 4B-K). Among the RNA-Seq data, SAP30, MT2A, 



sorting nexin 9 {SMX9), transthyretin (TTR), and matrix gla protein [MGPj, 
were up-regulated 126.47-fold, 17.77-fold, 5.75-fold, 2.5-fold, and 
2.1 -fold, respectively. RT-PCR also indicated an approximately 
four-fold increase in expression for SAP30, MT2A whereas TTR 
and MGP showed about 2 fold increase in their expression. 
However, RT-PCR results of SNX9 exhibited about 1.4-fold 
increase. Similarly, SCNIA, SYNC, RYRl, pkxin 1 [PLXMCl), and 
copine III [CPNE3) were down-regulated by 0.01 -fold, 0.02-fold, 
0.08-fold, 0.13-fold, and 0.17-fold (>4-fold), respectively. Howev- 
er, the RT-PCR data indicated an approximately two-fold 
decrease in their expression levels. 

Discussion 

The current study offers the first thorough insight into the 
transcriptome analysis of primary bovine MSCs with MYOGj^d 
using RNA-Seq technology. The number of total reads that map 
to the reference genome met the high quality criterion of the 
RNA-seq technology [64]. The most practical justification for 
reads not mapping uniquely to the reference genome could be due 
to the sequencing errors or polymorphisms, reads that come from 
repetitive sequences, and reads from exon-exon junctions [65] . 

Several genome wide high-throughput studies have been 
applied to investigate the functional role of various genes during 
myogenesis [4], [39], [66], [67]. Recendy, a microarray based 
study o{ MyoG has shown its role in mediating cell cycle exit in the 
absence oti p38oi and recognized an important function oi p38a. in 
cell fusion through the up-regulation of CD53 [68] . Another DNA 
microarray based study identified approximately 193 additional 
transcriptional regulators with varying expression levels during 
myogenesis [41]. DNA microarray has also been used to observe 
global changes in C2C12 cells transcriptome stimulated by 
exogenous myostatin (also known as GDF8) treatment, as well as 
to identify a network of genes involved in the inhibitory effects on 
dififerentiation [69]. In addition to microarray based studies, the 
RNA-Seq technique has been apphed using C2C12 mouse 
myoblast cell lines to detect 13,692 known transcripts and 3,724 
unannotated transcripts [29]. These sequencing or array-based 
methods have been shown to improve our understanding of 
myogenesis by revealing a broad range of target genes of myogenic 
transcription factors, novel myogenic factors and the characteris- 
tics of myoblasts and myotubes, which are difficult to identify by 
traditional approaches. 

However, almost aU of the aforementioned studies have used 
C2C12 mouse cell lines. We recently used primary bovine cells of 
high purity [70] to identify genes differentially expressed during 
differentiation and transdififerentiation of MSCs and differentia- 
tion of preadipocytes [37], [70]. MSCs are stem cells that reside 
between the sarcolemma and the basal lamina of adult skeletal 
muscle [71]. Since the serum derived from bovine species is an 
essential component of the in vitro cell culture system, there would 
be a great advantage of using bovine primary MSCs that closely 
mimic the in vivo situation during myogenesis [72]. Indeed, such 
studies might enable enhancement of muscle fiber characteristics, 
leading to improved meat quality. MyoG is specifically responsible 
for muscle fiber characteristics and closely associated with meat 
quality by affecting muscle development [73], [74]. 

Key processes altered during MyoGkd 

Genes involved in cell cycle regulation and DNA 
replication. MyoGj^j caused up-regulation of a large number 
of genes involved in functions related to cell proliferation, such as 
DNA replication, the cell cycle and mitosis (Table 3A & Table 
S3A). Figure 2B verified the computational results showing that 
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Table 3. Significantly enriched gene ontology terms detected by FAC in A) up-regulated genes, and B) down-regulated genes. 





S. No. 


GO Term (Fold enrichment) 


No. of Genes 


P-Value 


A) Up-regulated 


1. 


GO:0005654~nucleoplasm (1.92) 


22 


0.0047 




2. 


GO:0031 981 -nuclear lumen (1.54) 


29 


0.0189 




3. 


GO:0031974~membrane-enclosed lumen (1.45) 


35 


0.0203 




4. 


GO:0043233~organelle lumen (1.44) 


35 


0.0254 




5. 


GO:0070013~intracellular organelle lumen (1.43) 


33 


0.0306 




6. 


GO:0005615~extracellular space (1.91) 


17 


0.0157 




7. 


GO:0044421 -extracellular region part (1.68) 


21 


0.0227 




8. 


GO:0004222~metalloendopeptidase activity (3.85) 


5 


0.0404 




9. 


GO:0043933— macromolecular complex subunit organization (1.76) 


16 


0.037 




10. 


GO:0000280~nuclear division (2.84) 


8 


0.0225 




11. 


GO:0007067~Mitosis (2.84) 


8 


0.0225 




12. 


GO:0000087~M phase of mitotic cell cycle (2.79) 


8 


0.0245 




13. 


GO:0048285~organelle fission (2.73) 


8 


0.0272 




14. 


GO:0051 301 -cell division (2.39) 


9 


0.0349 




15. 


GO:0007049~cell cycle (1.71) 


17 


0.0384 




16. 


GO:0005819~Spindle (3.14) 


6 


0.0418 




17. 


GO:0005875~microtubule associated complex (3.77) 


5 


0.043 




18. 


GO:0000278~mitotic cell cycle (2.11) 


10 


0.0469 




19. 


GO:0007346~regulation of mitotic cell cycle (3.09) 


6 


0.0447 


B) Down-regulated 


1. 


GO:0004721 —phosphoprotein phosphatase activity (3.51) 


7 


0.0146 




2. 


GO:0006470~protein amino acid dephosphorylation (3.72) 


6 


0.0224 




3. 


GO:001 631 1 -dephosphorylation (3.21) 


6 


0.0387 




4. 


GO:0043292~contractile fiber (4.06) 


6 


0.0159 




5. 


GO:a030017~sarcomere (4.18) 


5 


0.0312 




6. 


GO:0030016~myofibril (3.69) 


5 


0.046 




7. 


GO:0044449~contractile fiber part (3.63) 


5 


0.0486 




8. 


GO:0000904~cell morphogenesis involved in differentiation (3.38) 


10 


0.0028 




9. 


GO:0a00902~cell morphogenesis (2.55) 


11 


0.0111 




10. 


GO:0032989~cellular component morphogenesis (2.29) 


11 


0.022 




11. 


GO:0031 175~neuron projection development (2.58) 


8 


0.0355 




12. 


GO:0048666~neuron development (2.19) 


9 


0.0529 




13. 


GO:0016049~cell growth (8.11) 


6 


0.0008 




14. 


GO:0040007~growth (3.61) 


8 


0.0067 




15. 


GO:0008361 -regulation of cell size (3.20) 


8 


0.0126 




16. 


GO:0032535~regulation of cellular component size (2.74) 


9 


0.017 




17. 


GO:0031090~organelle membrane (1.64) 


22 


0.0239 




18. 


GO:0031966~mitochondrial membrane (2.08) 


10 


0.0506 




19. 


GO:0005874~microtubule (2.69) 


9 


0.0186 




20. 


GO:0005856~cytoskeleton (1.48) 


25 


0.0445 




21. 


GO:0006796~ phosphate metabolic process (1.70) 


20 


0.0251 




22. 


GO:0006793~phosphorus metabolic process (1.70) 


20 


0.0251 




23. 


GO:0030424~axon (3.09) 


6 


0.0443 



doi:1 0.1 371 /journal.pone.0092447.t003 



CCNA2 expression increased by more than ten-fold in response to 
MyoGitd based on real time RT-PCR. Cell cycle related genes in 
these groups consist of several cell division homologue genes 
including cell division cycle 45 {CDC45), cell division cycle 20 {CDC20) 
and cell division cycle 6 {CDC6), each showing more than a four-fold 
change. Promoter studies in quiescent myoblasts have shown that 
MyoD activates the expression of CDC6 and minichromosome 



maintenance complex component 2 (A4CM2) genes, which prepare 
chromatin for DNA replication and as a result progression of the 
cell through the S-phase. Additionally, several other key cell 
division cycle-associated proteins including cell division cycle associated 
2 {CDCA2), cell division cycle associated 3 {CDCA3), cell division cycle 
associated 7 [CDCA7) and cell division cycle 8 [CDC8) were up- 
regulated in response to MyoGi;,!. MyoG plays a critical role in 



PLOS ONE I www.plosone.org 



8 



March 2014 | Volume 9 | Issue 3 | e92447 



RNA-seq Analysis of MyoG Knock-Down in Bovine IVlSCs 



A) 



B) 



C) 



D) 




on 0.0 



1(3 ^ o ^ ^ 



E) 



F) 



{: ^ ^ ^ !:: ^ i 

i M B B m I' l M m 

CC 0 ' ^ ^ a 0.0 J ^ ^ DC 0.0 ' — ^ 

MyoG„<, MyoGkd IVlyoG„d MyoG^j MyoG^,, MyoGk^ 

{p= 0.0016) (p= 0.0758) (p= 0.0184) 



c 3.0 
o 

m 2.5 

in 

<u 

£ 2.0 

2 1.5 

0) 

.2 1.0 
10 

■3; 0.5 
CC 

00 



I) 



1.2 
1.0 



£ 0.8 
a. 

S 0.6 



is 0.2 
o 

0.0 



MyoG„d MyoGkd 

(p= 0.0118) 



iVlyoG„ 



{p= 0.0016) 



G) 



MyoGk, 

(p= 0.0758) 
H) 



1 MGP g 1.6, SNX9 a I ''•2 I SCWM g 1-2 1 SVWC 

Ui i«l iliiii 

MyoG„d MyoGkd MyoG„<, MyoG,,,, MyoGw^ MyoG|<d MyoG„d MyoG|<d 

(p=0.1155) ,^ (p=0.0206) (p= 0.0377) {p= 0.0027) 

J) K) 

] RVRf a o 1 PLXNC1 g 1 CPNE3 

ii iii lit 



MyoG„ 



MyoGm 
(p=0.0002) 



MyoG„d MyoGkd 

(p=0.0004) 



Figure 4. Real time RT-PCR validation of muscle specific and differentially expressed genes on MyoGkd. A) Fold changes of Myf5, MyoD, 
MyoG, MYL2, and MYH3 genes determined by RNA sequencing were compared with real-time RT PCR results. B-K) RT-PCR validation of mRNA 
expression for five randomly selected genes confirmed the increased expression of MT2A SAP30, TTR, MGP, and SNX9 genes and decreased expression 
of SCNIA, SYNC, RYRl, PLXNCl, and CPNE3 genes by IVlyoGkd- MyoG„d represents control, respectively (mean ±S.D., n = 3). p-value indicates the 
statistical significance of the data and different letters (a and b) in graph show significant differences among groups. 
doi:10.1371/journal.pone.0092447.g004 



mediating terminal difTerentiation through cell cycle exit, and the 
activation of several cell cycle genes as a consequence of AfyoG 
down-regulation is well known [68]. Among the transcription 
factors, E2F transcription factor 1 {E2F1), which is known to play 
important roles in regulating cell proliferation [75], was up- 
regulated by about seven-fold. Previous studies have demonstrated 
that the expression of MjioG is strongly correlated with miRNA 
(miR-20a) expression, which in turn controls cell cycle exit by 
targeting E2F transcription factors [68], [76], [77], [78]. Two high 
mobility group box genes (HMG20B and HMGAT) that play a role 
in the regulation of DNA-dependent processes (transcription, 
replication, and DNA repair) involved in altering the conforma- 
tion of chromatin also belong to this group of up-regulated 
biological processes [79] . Moreover, expression of MCM proteins 
(putative replicative helicase) such as MCM5 is necessary for DNA 
replication [80], [81] and essential for DNA replication fork 
progression [82], [83]. 

Processes related to organelle lumen, nucleoplasm and 
cytosol. DAVID functional analysis identified 59 unique genes 
up-regulated by more than four-fold representing the GO terms 
organelle lumen, nucleoplasm and cytosol (Table 3A and Table 
S3A). These genes included transcription factor E2F1, endoplasmic 
reticulum resident protein 44 {ERP44), gamma-enolase [EJV02) and vascular 



endothelial growth faitor-D ( VEGE-D), which are involved in a wide 
variety of biological processes. The proteins belonging to the E2F 
family of transcription factors play a significant role in controlling 
cell proliferation. For example, E2E1 is considered the key target 
of pRB and is regulated by pRB throughout cell proliferation [84] . 
ERP44, a thioredoxin [TRXj family protein known to be involved in 
oxidative protein folding, directly regulates or inhibits the channel 
activity of inositol 1 ,4,5-trisphosphate receptors [IPsRs). ERP44 exclu- 
sively interacts with the L3V domain oi IPjRl, and this binding is 
dependent on pH, Ca * concentration, and redox state [85]. 
EJV02, a membrane protein, has been reported as a marker of 
neuro-muscular junctions [MMJs], whose expression decreases consid- 
erably during the initial stages of human embryonic muscle tissue 
development [86] . VEGE-D interacts with, and induces dimeriza- 
tion and tyrosine autophosphorylation of its endothelial cell- 
specific receptor, VEGER-2, which stimulates endothehal sprout- 
ing, proliferation, and survival, as well as vascular permeability. 
Similarly, binding of VEGF-D with VEGFR-3 stimulates related 
processes in lymphatic endothelial cells [87], [88], [89], [90]. 

Phosphorus metabolic process. FAC analysis identified 
phosphorous metabolic processes as an important biological 
process in MyoGj^d ceUs (Table 3B & Table S3B). Both the 
FAC and function annotation chart analysis detected about 21 
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unique genes that exhibited s4-fold decrease in their expression 
rates and were involved in phosphorous related processes. These 
genes included various phosphatases and kinases such as receptor- 
type tyrosine-protein phosphatase beta (PTPRB), serine / threonine-protein 
phosphatases PPl-beta catalytic submit {PPPICB), serine / threonine-protein 
kinase 40 {STK40) and cyclin-dependent kinase 13 {CDK13). One of the 
main mechanisms by which the signaling cascades control various 
stages of myogenesis is through protein kinases that direct cell 
behavior via the reversible process of phosphorylation [91]. 
Extensive studies have revealed that protein tyrosine phosphatases 
play a vital role in regulation of skeletal muscle myogenesis [59], 
[92], [93], while dephosphorylation of tyrosine residues is required 
for cell cycle exit during myogenesis [94]. The KEGG pathway 
analysis identified various pathways that lead to down-regulation 
of various protein phosphatases during at least one step of their 
representative pathway, such as the PI3K-Akt signaling pathway, 
MAPK signaling pathway, focal adhesion, TGF-beta signaling 
pathway and hippo signahng pathway (Figure 5). One of the 
down-regulated genes is protein phosphatase 1, catalytic subunit, 
beta isozyme {PPPICB), which encodes a serine/threonine-protein 
phosphatase PPl-beta catalytic subunit, an important enzyme 
responsible for protein phosphorylation and regulation of many 
physiological processes [95] . Pathway analysis also illustrated that 
PPPICB is involved in many important pathways related to 
myogenesis such as focal adhesion, the Hippo signaling pathway 
and regulation of the actin cytoskeleton (Figure S3A-E). 

Cytoskeleton and cell morphogenesis. DAVID FAC 
indicated that the GO terms cytoskeleton and cell morphogenesis 
involved in differentiation were down-regulated in response to 
MyoGkd. The analyses identified about 39 unique genes involved 
in these processes that showed S 4-fold reduction in their 
expression rates. The location of most of the genes in this category 
is either cytoplasm or labeled as secreted in the Uniprot database. 
The genes under these biological processes are involved in a broad 
range of functions su[:h as signaling pathways, transport, 
differentiation, etc. Some of these genes include disabled homolog 2 
(DAB2), microtubule-associated protein 2 [MAP2), synaptopodin-2 or 
myopodin (SXNP02), and moesin (MSN). Among these genes, DAB2 
(expressed in various tissues), which is detected at an early 
myogenic differentiation state [67], has lost or reduced expression 
in hyperproliferative cells [96] . Another gene in this group that is 
significandy down-regulated is a member of the tissue inhibitors of 
matrix metalloproteinases (TIMP) family, TIMP3, or metalloprotei- 
nase inhibitor 3. TIMP3 complexes with MMPs and is the only TIMP 
capable of inhibiting membrane bound MMP, transmembrane 
MMP and sheddases such as TNF-ot converting enz)rme (TACE), 
which is also known as disintegrin and metalloproteinase (ADAM- 
17) [97], [98]. Conversely, all MMPs detected in this study were 
highly up-regulated. Pathway analysis also revealed that one of the 
ERM proteins known to regulate cross-linking of the plasma 
membrane and actin cytoskeleton, MSN, was down-regulated [99], 
[100], [101]. 

Role of glycogene expression in myogenesis 

Skeletal muscle development consists of a well controlled and 
regulated progression of various cellular processes, including cell 
proliferation, migration, and differentiation [102]. Until recently, 
the role of glycoproteins in myogenesis did not receive a great deal 
of attention from the scientific community [5], [102]. However, 
many independent studies have recendy focused on the numerous 
roles of glycoconjugates during myogenesis [102]. The results of 
these studies have indicated that the expression of MyoG is pardy 
regulated by the reduced glycosylation-dependent recruitment of 
MeJ2D to MyoG promoter, suggesting negative regulatory mech- 



anisms of skeletal muscle development by O-GlcNAc glycosylation 
[103]. 

Different processes related to the formation and maintenance of 
skeletal muscles are characterized by the expression of a wide 
variety of molecules that strongly alter biological events, such as 
muscle development, differentiation and regeneration. Among the 
different types of macromolecules participating in myogenesis, 
interest in glycoproteins has been gaining remarkable attention in 
recent years; however, there are still several unanswered questions 
regarding their roles during skeletal muscle development [102]. 
Similar to other eukaryotic cells, the plasma membrane and ECM 
of myoblasts are rich in glycoproteins and glycolipids [5]. 
Inhibition of some ECM proteoglycans [syndecans) has shown to 
stop the progression of myoblast proliferation and fusion, 
regardless of the expression of MRFs [104], [105]. Similarly, 
interrupting N-glycan synthesis affects myoblast fusion [106]. 
Glycolipids also play key roles in cell differentiation and muscle 
development [107], [108]. 

Role of channels in myogenesis 

The initial phases of myogenesis are marked by the develop- 
ment of excitability and contractile properties by skeletal muscle 
cells [109]. Voltage dependent sodium channels comprise one of 
the key types of proteins that play a pivotal role in propagating 
action potential in nerves and muscle [110], [111], [112]. SCNIA 
consists of four homologous domains [113], and its activity is 
regulated by the interaction with fibrobhst growth factor 13 {FGF13) 
[114]. Many mutational studies have recognized a variety of 
tainted channel properties that include changes in the voltage 
dependence of activation and inactivation, speedy recuperation 
from inactivation, improved constant current and loss of channel 
function [113]. Similarly, RTRl, a calcium channel that plays an 
important role in excitation-contraction coupling in skeletal 
muscles, has shown increased expression levels during the early 
stages of myogenesis along with dihydropyridine receptors [DHPRs) 
[115]. STNC, which belongs to the intermediate filament (IF) family of 
proteins [116], is greatiy expressed in skeletal and cardiac muscles 
and may play an important role in maintaining contractile 
properties [117]. STNC is known to interact with another member 
of the IF family, desmin, and may play a significant role in 
protecting muscle cells from mechanical stress [118], [119] in a 
fashion similar to that of other members of the IF family [120]. 

Together, these data offer extensive and deeper insight into the 
transcriptional regulation of myogenesis by MyoG and provide a 
rich scope for designing future experiments to elucidate the 
pathways involved in skeletal muscle differentiation. Furthermore, 
to improve muscle growth as well as quality and quantity of meat, 
it is essential to recognize how MyoG influences the differentiation 
of MSCs. 

Conclusion 

In summary, our transcriptome analysis of primary bovine cells 
using RNA-Seq offers important insight into the transcriptional 
regulation of gene expression in down-regulated MyoG muscle 
cells. In addition to the identification of new genes in skeletal 
muscle development, our bioinformatics analysis suggested a role 
of phosphorous metabolic processes, proteins with channeling 
function, and the involvement of a significant number of 
glycogenes in my()g(;nesis. Further investigation of the genes 
identified in this study wiU facilitate our understanding and help 
explain the mechanism responsible for increasing skeletal muscle 
mass. 
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Figure 5. KEGG pathway map analysis of differentially expressed genes (&4-fold change). Among the various different pathways 
reported by KEGG analysis, two representative pathways A) PI3K-AKT and B) MARK are shown. For each rectangular shape, red and green borders 
indicate up-regulated and down-regulated genes, while the yellow star indicates a role of phosphoprotein. Each KEGG pathway analysis figure 
depicts the role or involvement of a specific gene or a family of related genes at a particular location in the pathway. For instance, the FGF gene in 
KEGG MARK pathway map represents the other FGF members (such as FGF10, FGF1 1, FGF12, etc) as well that may or may not be affected. Here, FGF 
represents all the members of fibroblast growth factor family. For FGF family members, the gene that exhibited >4 fold increase includes FGF11 
whereas FGF10, FGF16, and FGF18 showed >1 fold increase in their expression rates. Conversely, the other FGF members that showed <4 fold 
decreased expression include FGF12 and FGF14. Similarly, FGFR represents the other members (FGFR1, FGFR2, etc) of this family of receptors as well. 
In case of FGFR family, FGFR2 showed >4 fold decrease in its expression, and FGFR4 exhibited >1 fold decrease in its expression rate. However, 
FGFR1 and FGFR3 showed about 1 fold increase in their expressions. 
doi:l 0.1 371 /journal.pone.0092447.g005 



Supporting Information 

Figure SI Pax7 expression in MSCs. Cellular localization of 
Pax7 on MSCs at Day 10 by immunocytochemistry. A) Cell 
picture B) DAPI-stained nuclei. C) Pax7 antibody stained cells. 
(TIF) 

Figure S2 Multiple sequence alignment of LYLl and 
other bHLH genes. Multiple sequence alignment of LYLI 
(UNIPROT ID: E1BAR3) protein and other bHLH proteins 

(MyoD, MyoG, Myf5, and herculin) which are involved in skeletal 
muscle development. The region in the box indicates high 
sequent e similarity in the bHLH region of these muscle specilic 
proteins. 
(TIF) 

Figure S3 Other muscle specific pathways affected by 

MYOGkd. A) WNT signaling pathway, B) Focal adhesion, C) 
TGF-beta signaling pathway, D) Hippo signaling pathway, E) 
PPAR signaling pathway. 
(PDF) 

Table SI Primer information. 

PCLSX) 
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